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INTRODUCTION 

The  presence  of  pressure  fluctuations  on  a  turbomachine  blade  row  is 
known  to  cause  undesirable  effects,  such  as  vibration,  radiated  noise,  and 
performance  degradation.  It  is  important,  therefore,  to  provide  a  means  by 
which  a  turbomachine  designer  can  predict  the  unsteady  response  of  a  blade 
row  and  select  proper  design  parameters  to  minimize  these  effects. 

The  present  analysis  represents  an  alternative  method,  which  is  derived 
from  the  previous  theory  developed  by  Henderson  [1],  to  determine  the  unsteady 
lift  in  a  turbomachine  blade  row.  The  flow  is  assumed  two-dimensional,  incom¬ 
pressible,  and  inviscid.  The  vortex  representation  of  the  unsteady  cascade 
flow  is  shown  in  Figure  2  in  Reference  1.  An  explicit  solution  has  been 
derived  for  the  unsteady  pressure  difference  on  the  surface  of  a  turbomachine 
blade  row  operating  in  a  spatially  varying  disturbance  flow  field.  Figure  1. 
Specifically,  this  solution  is  expressed  in  terms  of  the  cascade  geometrical 

parameters — space-chord  ratio,  s/c,  stagger  angle,  £,  blade  camber,  y  ,  mean 

m 

angle  of  incidence,  a^,  and  the  disturbance  flow  characteristics,  represented 
by  the  reduced  frequencies,  U)  and  X. 

One  of  the  advantages  of  employing  the  present  analysis  is  that  the 
parameters  of  unsteady  lift  and  pitching  moment  can  be  easily  obtained  by 
direct  numerical  integration  of  the  unsteady  pressure  difference  across  the 
airfoil  chord. 

A  program  has  been  written  to  formulate  the  unsteady  pressure  distribution 
on  a  single  blade  of  a  blade  row  as  a  function  of  the  design  parameters  and  it 
performs  an  integration  of  the  pressures  to  obtain  the  unsteady  lift  and 
pitching  moment  using  a  quadrature  method. 

The  significance  of  the  present  program  can  be  summarized  as  follows: 

1.  Because  the  effect  of  blade-to-blade  interaction  is  considered,  the 
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resonance  phenomena  which  occurs  when  the  blade  spacing  equals  the 
disturbance  wavelength  can  be  predicted. 

2.  The  effect  of  the  blade  camber  in  a  cascade  of  airfoils  is  included. 

3.  The  effect  of  the  mean  incidence  angle  of  the  airfoils  is  included. 

4.  The  program  accounts  for  either  a  nonconvected  or  convected  disturbance. 

5.  The  program  gives  a  complete  description  of  the  unsteady  response,  i.e., 
the  unsteady  pressure  distribution,  lift,  moment,  and  center -of -pressure 
location. 

DEFINITION  OF  THE  PROBLEM 

The  problem  considered  is  the  unsteady  pressure  distribution  on  the  blades 

of  a  cascade  when  the  blade  row  experiences  a  spatial  variation  in  the  inflow; 

for  example,  the  nonuniform  flow  caused  by  the  wakes  of  an  upstream  blade  row. 

The  flow  field  is  assumed  to  be  two-dimensional,  inviscid,  and  incompressible, 

thus  representing  the  development  of  a  cylindrical  surface  as  shown  in  Figure  1. 

This  represents  the  passage  of  a  rotor  through  the  wakes  of  an  upstream  stator 

blade  with  a  swirling  mean  flow.  These  wakes  have  a  maximum  velocity  deficit  of 

2w^  and  are  transported  over  the  rotor  by  the  velocity,  W^,  the  circumferential- 

mean  velocity  relative  to  the  rotor  blades  with  the  wake  present.  The  velocity 

deficit,  w , ,  represents  the  fluctuation  about  the  mean  velocity,  W  . 

d  m 

The  description  of  the  wake  deficit  shown  in  Figure  1  can  be  accomplished 
by  using  Fourier  series  representation.  From  this  representation,  the  contri¬ 
bution  of  each  harmonic  of  the  velocity  variation  to  the  unsteady  response  of 
the  blade  can  be  determined.  The  resulting  components  of  the  disturbance 
velocity,  w^,  normal  and  parallel  to  the  chord  can  be  expressed  from  the 
geometry  of  Figure  1  as, 

v,  *  -  w,  sin  4> 
d  d 


u ,  ■  -  w.  cos 
d  d 


(1) 
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The  consideration  of  the  disturbance  as  components  parallel  and 
normal  to  the  chord  facilitates  the  solution  of  this  problem  and  it  is 
possible  to  add  the  respective  solutions  since  the  model  employed  is 
linear . 

FORMULATION  OF  THE  UNSTEADY  PRESSURE  DISTRIBUTION 

In  order  to  generate  values  of  the  coefficients  of  unsteady  pressure 
difference,  AC^,  unsteady  lift,  and  unsteady  moment,  C^,  it 
is  necessary  to  consider  the  contribution  of  blade  geometry  to  this  unsteady 
response.  Starting  from  Equation  (A20)  in  Reference  1,  and  substituting  the 
boundary  conditions,  Equation  (14),  into  this  equation,  the  function  of 
unsteady  pressure  distribution  can  be  determined  and  is  expressed  in  an 
explicit  form  as: 


(2) 
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The  quantities  A  and  y  (a  )  are  defined  as 

s 

2ueiU,{y^d[J0(X)  -  J2(X)  -  2J1(X)]  +  <?d[JQ(X)  -  iJ^X]} 


,iUr„(2), 
1 


ia)ireiW[H1v‘-' (tu)  +  1hJ2)(u))]  +  M_ 


-  1  oo 

-  1  oo 

M  = 
c 

l  +  l 

-CO  1 

inTr„  ,  _  .  io) 

e  [C^  +  C2  -  2]  -  lore 

l  +  l 
-00  1 

einX[D1  +  D2] 


(3) 


and 


Y0  (o+> 
s 


=  2W  / 


■'  l-o 


m 


/  1  +  o' 


C(i+o+> 


m 


y++a  ) 

'm  m 


-1« 

l+l 
-°°  1 

(E^Ej) 

2  - 

-  1  00 
l+l 

_oo  1 

(C1+€2-  2) 

• 

(4) 


for  a  cascade  which  experiences  a  nonconvected  disturbance  of  reduced 
frequency  A.  The  flow  angle,  $>,  is  related  to  the  stagger  angle,  £,  and 
the  mean  swirl  angle,  8,  entering  the  cascade,  as  <J>  =  it  -  g  - 

Equations  (2) ,  (3) ,  and  (4)  contain  several  quantities  expressing  the 
contribution  of  the  blades  in  a  cascade  which  are  adjacent  to  the  reference 
blade.  Each  of  these  quantities  are  infinite  summations  and  are  listed  in 
Table  I. 


NUMERICAL  TECHNIQUES 

Several  numerical  methods  were  employed  to  give  a  solution  to  the 
problem  described  above.  They  are  discussed  in  the  following. 


Averaging  Procedure 

The  evaluation  of  the  infinite  summations  having  the  form 


A 


in. the  present  analysis,  as  listed  in  Table  I,  uses  the  convergence  criterion 
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called  "Averaging  Procedure."  This  method  was  previously  used  by  Henderson  [1] 
In  his  computer  program  to  calculate  some  of  these  unsteady  cascade  terms. 

The  mathematical  expression  of  the  criterion  is: 


D1FF  = 


N+l 

-  1  N 

N 

- 1  N 

l 

z  +  z 

A  I 

z  +  z 

3 

-N  1 

3 

-  N  1 

N  -  2  N  -  3 


<  e 


(5) 


The  number  of  pairs  of  arguments,  N,  required  for  an  accurate  estimation  of 
these  infinite  cascade  summations  can  be  determined.  It  was  observed  that 
as  the  value  of  e  decreases,  the  number  of  argument  pairs  which  are  summed 
increases  and  the  summation  approaches  a  finite  value  in  an  asymptotic  manner. 

An  example  is  illustrated  in  Figure  2.  The  validity  of  the  method  is  there¬ 
fore  assured  by  this  converging  behavior. 

Gauss-Legendre  Quadrature  Method  ' 

When  evaluating  the  unsteady  lift  and  pitching  moment  using  the  calculated 
unsteady  pressure  difference  across  the  airfoil  chord,  a  difficulty  arises  when 
performing  the  integration  from  the  leading  to  trailing  edge  due  to  the  exis¬ 
tence  of  a  singularity  at  the  leading  edge.  This  difficulty  was  overcome  by 
the  use  of  the  Gauss-Legendre  (G-L)  method.  The  Gauss-Legendre  quadrature 
employed  has  a  behavior  which  is  most  helpful  with  the  present  problem,  and,  in  addi¬ 
tion,  the  very  high  order  interpolating  polynomial  which  is  implicit  in  Gauss 
quadrature  can  better  approximate  the  value  of  integration  near  the  singu¬ 
larity  [2].  The  region  of  integration  is  further  divided  into  two  parts  such 
that  one  of  the  integrals  is  carried  out  over  a  small  region  close  to  the 
singularity.  For  example, 
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1 

SL  -ST  I  “,(«+>"+ 

-1 

-.98  1 

=  |  ACp(o+)do+  +  |  ACp(a+)da+ 

-1  -.98 

The  theory  by  Sears  for  the  unsteady  response  of  an  Isolated  flat- 
plate  airfoil  [3]  wasused  t0  check  the  accuracy  of  this  integration  scheme. 
Excellent  agreement  was  found  between  the  results  using  the  G-L  method 
and  those  obtained  by  the  exact  solution  developed  by  Sears,  Figure  3. 

COMPUTER  PROGRAM 

The  expression  developed  for  unsteady  pressure  distribution  is  intro¬ 
duced  above.  The  Equations  (2),  (3),  ?nd  (4)  have  been  programmed  using 
complex  Fortran  IV  language.  This  program  can  be  run  on  the  P.  S.  U. 

360  system.  Note  that  memory  size  S  *  280K  must  be  used. 

The  computer  program.  Appendix  A,  consists  of  a  main  program  and 
several  subroutines.  The  arrangement  of  this  program  and  its  subroutines 
is  as  follows: 

MAIN  PROGRAM 
Subroutines  UPDBC 

CCSUM 
RSUM1 
RSUM2 
BESJ 
BESY 
ZED 


ATANF 


JP. 7T" 
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The  main  program  also  performs  the  integration  required  to  find  unsteady 

i 

lift  and  moment.  Part  of  the  major  effort  was  to  generate  the  subroutines 
permitting  the  evaluation  of  the  infinite  summations  occurring  in  the  t 

analysis.  The  subroutines  are  listed  and  described  below: 


SUBROUTINES 

INPUT 

OUTPUT 

DESCRIPTION 

UPDBC 

s/c,C,y  , 
a  , 

m 

us,  A 

a+ 

AC 

s 

AC 

s 

The  calculation  of  unsteady  pressure  difference 
at  certain  locations  along  the  airfoil  chord  and 
its  components  contributed  by  chordwise  and 
transverse  disturbances  u<j  and  v^,  respectively. 

CCSUM 

s/c, 

HI, 

a+, 

e 

-  1  N 

M 

-N  1 

Evaluation  of  the  infinite  cascade  summations 
having  the  form  in  Table  I.  Types  of  the 
summations  are  controlled  by  the  flag,  1, 

I  -  1,  2,  3,  ...  7.  ; 

] 

RSUM1 

aj,  ns/c, 

K 

Dl+D2 

i 

Evaluation  of  the  integrals  representing  the 
contribution  due  to  the  wakes  shed  by  nth 
blades. 

RSUM2 

03,  ns/c, 

Bl  +  B2 

Evaluation  of  the  integrals  representing  the 
contribution  due  to  the  wakes  shed  by  nth 
blades. 

BESJ,  BESY 

03, m 

Jm’Ym 

Determines  Bessel  functions  of  the  first  and 
second  kind,  order  m  and  argument  03. 

ZED 

A,o+ 

Z 

A  function  resulting  from  the  convected 
disturbance  effect. 

ATANF 

R,  l 

0 

Conversion  of  the  vectors  into  polar  coordinates. 

Since  the  present  analysis  determines  the  effects  of  blade-to-blade 


interactions,  the  program  is  also  capable  of  solving  the  isolated  airfoil 
case.  As  the  value  s/c  increases,  values  of  the  cascade  terms  listed  in 
Table  I  are  found  to  decrease.  If  s/c  is  large  enough,  their  contributions 
would  become  negligibly  small.  To  save  computational  time,  the  program 
automatically  treats  these  unsteady  cascade  terms  as  zero  when  the  space- 
chord  ratio,  s/c,  is  set  to  a  number  larger  than  100.0. 
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For  a  factor  of  e  -  0.001,  It  takes  about  ten  seconds  to  perform  the 
calculations  for  each  value  of  w.  Care  must  be  taken  for  the  cases  where 

values  of  to  are  in  the  neighborhood  of  resonance  points  or  are  extremely 
small.  In  these  cases,  the  computational  time  requited  to  obtain  a  solution 
becomes  significantly  long. 


PREPARATION  OF  THE  INPUT  DATA  CARDS  AND  EXAMPLES 

The  computer  program  must  be  arranged  in  the  order  described  in  the 
previous  section.  To  activate  the  unsteady  pressure  program,  the  input 
information  must  be  arranged  on  four  input  cards  in  the  following  fashion: 


COMPUTER 

SYMBOL 

SYMBOL  IN 
ANALYSIS 

DESCRIPTION 

REMARKS 

CARD  1  (3F10.4) 

XC 

+ 

X 

c 

Dimensionless  location  of 
concentrated  vortices  on 
neighboring  blades. 

x*  =  -0.50  in  the  analysis 

ERR 

e 

Parameter  controlling  the 
accuracy  of  calculations 
of  cascade  summations. 

e  “  0.001  recommended. 

BETA 

e 

The  flow  exit  angle  follow¬ 
ing  front  stage  cascade. 

8=0°,  if  inf low  has  no 
swirl,  (degrees) 

CARD  2  (3110) 

KRUN  Number  of  runs  to  calculate  using  different  combinations 

of  parameters  described  on  Card  3. 

IRUN  Number  of  runs  to  calculate  using  different  values  of 

reduced  frequencies,  w  and  X,  which  are  described  on 
Card  4. 

JRUN  If  JRUN  *=  0,  print  out  C^  and  CM  only; 

«  1,  print  out  complete  results  including 

ACp  as  a  function  of  x*,  CL  and  CM< 
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CAPD  3 

(4F10.4) 

SOC 

s/c 

Space  -to-chord  ratio  of 
cascade. 

If  s/c  larger  than  100.0,  1 

the  program  goes  to 
isolated  airfoil  1 

calculation. 

EPS 

£ 

Stagger  angle  of  cascade. 

(degrees) 

YM 

+ 

y 

■'max 

Dimensionless  maximum  camber; 
ratio  of  maximum  camber  to 
the  blade  half -chord. 

AM 

a 

m 

Mean  incidence  angle  of  the 
flow  relative  to  cascade. 

(degrees) 

CARD  4 

(2F10.4) 

OMEGA 

U) 

Reduced  frequency  based  on  the 

mean  flow  velocity  (convected) . 

OMEGB 

A 

Reduced  frequency  based  on  the  transport  velocity 
of  disturbance  over  the  airfoil  (nonconvected) . 

To  illustrate  the  use  of  the  program  listed  in  Appendix  A,  the  following  | 

examples  are  presented.  1 

Example  No.  1 

■ 

To  preform  calculations  for  a  given  cascade  at 

several  values  of  reduced  frequency 

e  -  o° 

,  s/c  = 

2.209,  £  -  *5°,  ym 

=  0.05,  a  =  1.2°  ; 

in 

(i) 

a)  =  1.0, 

A  =  1.0; 

(ii) 

ai  =  2.0, 

A  =  2.0; 

(iii) 

a)  =  3.0, 

A  =  3.0. 

The  input  data  cards  will  be  arranged  as  shown  in  Figure  4. 


Example  No .  2 


To  perform  calculations  with  different  values  of  s/c,  £,  y  ,  or  a^,  for 
g  »  30°  and  several  values  of  reduced  frequency: 


(i) 

s/c  =  2.029, 

£  -  45°, 

ym 

-  0.05, 

a 

m 

»  1.20°; 

(ii) 

s/c  =  2.029, 

£  -  55°, 

ym 

=  0.05, 

a 

m 

*  3.50°;  - 

(iii) 

s/c  **  1.135, 

£  =  35°, 

ym 

**  0.0  , 

a 

m 

*  0.0°;  (flat  plate airfoiH 

with 

nonconvected  reduced 

frequencies  of: 

(i) 

a)  -  1.0,  A  - 

0.8; 

(ii) 

<*)  ■  2.0,  A  ■ 

1.7. 

Note, 

there  will  be  3  X  2 

different 

cases 

calculated . 

The 

input  data  cards 

will  be  arranged  as  shown  in  Figure  5. 
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Example  No.  3 

To  perform  calculation  for  an  isolated  airfoil*  i.e.,  neglecting  the  blade- 
to-blade  interactions,  while  still  including  the  effects  of  blade  camber  and  angle 
of  incidence,  set  s/c  to  a  number  larger  than  100.0.  Consider  the  case 

$  =•  30“ ,  s/c  -  -  45°,  y  -  0.08,  a  -  -1.2° 

m  ra 

(i)  a)  ■  0.50,  X  =  0.40; 

(ii)  w  -  1.00,  X  =  0.80. 

The  input  data  cards  for  this  example  are  listed  in  Figure  6. 
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Difference  using  the  Gauss-Legendre  Method 
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//DATA. INPUT 

DD  * 

-0.5000 

0.0010 

0.0000 

1 

3 

1 

2.0290 

45.0000 

0.0500  1.2000 

1.0000 

1.0000 

2.0000 

2.0000 

3.0000 

3.0000 

Figure  4 . 

Input  cards  for  Example  No. 

//DATA. INPUT 

DD  * 

-0.5000 

0.0010 

30.0000 

3 

2 

0 

2.0290 

45.0000 

0.0500 

1.2000 

1.0000 

0.8000 

2 . 0000 

1.7000 

2.0290 

55.0000 

0.0500 

3.5000 

1 . 0000 

0.8000 

2.0000 

1.7000 

1.1350 

35.0000 

0.0000 

0.0000 

1.0000 

0.8000 

2.0000 

1.7000 

Figure  5.  Input  cards  for  Example  No.  2 
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//DATA. INPUT  DD  * 


-0.5000 

0.0010 

30.0000 

1 

2 

1 

100.0000 

45.0000 

0.0800 

0.5000 

0.4000 

1.0000 

0.8000 

-1.2000 


Figure  6.  Input  cards  for  Example  No.  3. 
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C  SOURCE  LANGUAGE:  FORTRAN  IV 


* ***************************************************************** 

UNSTEADY  PRESSURE  PROGRAM 
WRITTEN  BY  I.  C.  SHEN  9/9/79 

CONTENTS 

MAIN  PROGRAM 
SUBROUTINES  UPDBC 

CCSUM 
RSUM1 
RSUM2 
BESJ 
BESY 

FUNCTIONS  ZED 

ATANF 

***************************** **********************************1*** 


/  EXEC  FWCG.PARM-' NOSOURCE, ISNS' 

/SYS IN  DD  * 

IMPLICIT  REAL*8(A-H,0-Z,$) 

REAL* 8  NSC.JKl  ,JK2,JK3,JK4 

COMPLEX*! 6  EAVA,EAVB,AVL , AVM,AVP , AVJ , AVK, AVER, 

+  ZED  Z  TPD  PDU  PDV 

+  HAn£o,’haM1,m£,NcJdELU,DELV,DEL,RI,EXW,EXL, 

+  CL , CLU , CLV , CPU , CLWH , CM , CMU , CMV 

DIMENSION  XP(24),TPD(24).PDU(24),PDV(24),XI(12),W(12) 

COMMON  / AVIV/  TAU , SINXI .COSXI 
COMMON  /RSUM/  JKI  ,JK2 ,JR3,  JK4 

COMMON  /MAIN/  EAVA.EAVB,AVL,AVM,AVP,AVJlAVK..Z,XC,ERR 
COMMON  /DELTA/  RJAO.RJAI . RJA2 . RYO , RY 1 , RJB0.RJB1 ,RJB2,HANK0,HANK.l , 
I  MC , NC , DEL , DELU , DELV , RI , EXW ,PHI , SINPH , COSPH 


WEIGHTING  COEFFICIENTS  FOR  THE  GAUSS-LENGENDRE  METHOD 
M=24  IN  EACH  INTERVAL  OF  INTEGRATION. 


DATA  PI/3. 1415927D0/ 

DATA  XI/Q.0640568929D0,  0.191 1188675D0 ,  0.3150426797D0, 

+  0.4337935076D0,  0.5454214714D0,  0.6480936519D0, 

+  0.7401241916D0  0.8200019860D0 ,  0.8864155270D0, 

+  0.9382745520D0,  0.9747285560D0,  0.9951872200D0/ , 

+  W/0. 1 27938 19 53D0 ,  0.1258374563D0,  0.1216704729D0, 

+  0.1 155056682D0,  0. 1O744427O1D0 ,  0.0976186521D0, 

+  0.0861901615D0  0.0733464814D0  0.0592985849DO, 

+  0.0442774388D0’  0.0285313886D0,  0.0123412298D0/ 

C 

C  READ  DATA 

c  *  ***************************************************** ************** 


C  *  SOC - SPACE-CHORD  RATIO  * 

C  *  EPS - STAGGER  ANGLE  (DEGREES)  * 

C  *  YM  - RATIO  OF  MAXIMUM  CAMBER  TO  HALF-CHORD  * 

C  *  AM  - MEAN  FLOW  INCIDENCE  ANGLE  (DEGREES)  * 

C  *  KRUN — NUMBER  OF  DIFFERENT  CASCADE  GEOMETRIES  * 

C  *  TO  BE  CONSIDERED.  * 

C  *  IRUN— NUMBER  OF  REDUCED  FREQUENCIES  TO  BECONSIDERED.  * 

C  *  JRUN  -  0  PRINT  OUT  ONLY  THE  UNSTEADY  LIFT  * 

C  *  AND  PITCHING  MOMENT.  * 

C  *  -  1  PRINT  OUT  ALL  THE  RESULTS.  * 


Q  ******************************************************************** 

15250 
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READ(5 ,14)  XC, ERR, BETA 

14  FORMAT(3Fl0.4) 

READ(5 ,15)  KRUN.IRUN.JRUN 

15  FORMAT (3110) 

DO  20  K-l.KRUN 

READ (  5 , 1 0  J  SOC ,EPS,YM,AM 
10  FORMAT(4F10.4) 

WRITE '  ‘ 

600  FORMA 
14X 


ITE(6 ,600)  S0C.XC.ERR.EPS, BETA, YM, AM 

,8X  'ERR-'  E12.4// 

14X, 'STAGGER-' , F10.4, '(DEG)' ,8X, BETA-' .F10.4 , ' (DEG) '// 
14X. 'YM-' .F10.4 ,8X, 'AM-' ,F10.4. ' (DEG) ' /)/) 


A********* 'j/ 


ill—  .1  IVlTlUA.  m-l-  .IIUIT, 

IF(SOC.GE.IOO.ODO)  WRITE(6,103) 

GO  TO  104 

103  FORMAT(////8X. ***********  ISOLATED  AIRFOIL  CASE 
120X, 'S/C  -->  00') 

104  PHI-180. DO-EPS-BETA 
COSPH-DCOS(PHI/57 .29547D0 ) 

SINPH»DSIN(PHI/57.29547DO) 

AM-AM/ 57.29547D0 
SINXI=DSIN(EPS/57.29547D0 ) 

COSXI-DCOS (EPS/5 7 .2954 7D0) 

DO  25  I=*1,IRUN 

READ (5, 40)  OMEGA, OMEGB 
40  F0RMAT(2F10.4) 

OMGAB-O MEGA/ OMEGB 

- PARAMETER  TESTING  THE  PRESENCE  OF  NON-CONVECTIVE  EFFECT. 

C-YM*YM+AM*AM 

- PARAMETER  TESTING  THE  PRESENCE  OF  STEADY  VORTICITY 

[EQN  IN  SHEN'S  THESIS) 

TAU— 2 . *SOC*OMEG A*DCOS ( BETA/5  7 . 2  954  7  DO ) /S INPH 
WRITE(6 ,650)  OMEGA, OMEGB, TAU 

650  FORMAT ( / // 2X# * ****^*****w***************************************** 

I******************** 

1 '  /  / 1 5X . 'OMEGA-' ,F10.4,8X, 'LAMDA-' .F10.4/15X, 'TAU-', F10.4//2X, 

1 ' ***************************************************************** 
1******'//) 

******************************************************************** 

DETERMINE  UNSTEADY  CIRCULATION  COEFFICIENT  (DEL)  BY 
EVALUATING  CASCADE  SUMMATIONS  THAT  ARE  NOT  FUNCTIONS 
OF  XP. 

******************************************************************** 


IF(SOC-100. )  102,101,101 

101  EAVA=DCMPLX(O.ODO,O.ODO) 

EAVB=DCMPLX(0 .0D0,0.0D0 ) 

AVL=DCMPLX(0.0D0,0.0D0) 

AVM-DCMPLX(0.0D0,0.0D0) 

GO  TO  120 

102  CALL  CCSUM( 1 ,OMEGA,0 .0D0 ,XC ,SOC,EPS,ERR,EAVA, IAVA) 
EAVB=DCMPLX(O.DO,0.DO) 

IF(C.EQ.O.DO)  GO  TO  107 

CALL  CCSUM(2, OMEGA, 0.0D0.XC, SOC, EPS, ERR, EAVB,IAVB) 
107  CALL  CCSUM(3, OMEGA, 0.0D0,XC, SOC, EPS, ERR, AVL.IAVL) 
CALL  CCS  UM(4 .OMEGA. 0 . 0D0 , XC , SOC , EPS , ERR , AVM , IAVM) 
120  RI=DCMPLX(0.0D0,1.0D0) 

EXW-CDEXP ( RI*0MEGA ) 

EXL=CDEXP(-RI*0MEGB) 

CALL  BESJ  (OMEGB, 0 ,RJB0 , 1 .OD-5, IEO 
CALL  BESJ  (OMEGB, 1RJB1.1. OD-5, IE1 
CALL  BESJ  (OMEGB, 2, RJB2,1. OD-5, IE2 
CALL  BESJ  (OMEGA, 0,RJA0,1. OD-5, IAO 
CALL  BESJ  ( OMEGA, 1,RJA1,1. OD-5, IA1 
CALL  BESJ  (OMEGA, 2, RJA2,l. OD-5, IA2 
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CALL  BESJ  (OMEGA, 2 ,RJA2 ,1 .0D-5, IA2) 

CALL  BESY  ( OMEGA, 0, RYO, IYO) 

CALL  BESY  (OMEGA, 1 ,RY1 , IY1) 

HANK0=DCMPLX(RJA0,-RY0) 

HANK1=DCMPLX(RJA1 ,-RYl ) 

MC-EAVA-RI*OMEGA*EXW*AVL 
NC=RI*OMEGA*PI*EXW*(HANK1+RI*HANKO)+MC 
DELU=*-2 .  *P  I*EXW*  (-YM*COSPH*  (RJB0-RJB2-RI*2 .  *RJB1 )  )  /NC 
DELV=-2.*PI*EXW*(-SINPH*(RJB0-RI*RJB1))/NC 
DEL=DELU+DELV 
C 

Q  ****************  *******************  ***************  *****  ******** 

C  SEPARATE  THE  INTERVAL  OF  INTEGRATION  FROM  XP+-  [-1,  1]  INTO 
C  XP+  =  [-1,  -0.98]  AND  [-0.98.  1]  WITH  24  UNEVEN  STEPS  IN  EACH 

C  REGION  TO  GIVE  AN  ACCURATE  ESTIMATION  OF  THE  IMPROPER  INTEGRAL. 

C  *************************************************************** 

c 

M-0 

11  DO  33  J-1.12 

XP ( J ) — . 99D0- . 0 1DQ*XI ( 1 3-J ) 

33  XP(J+12)=-.99D0+.01D0*XI(J) 

GO  TO  222 

12  DO  22  J-1,12 
XP(J)=.01D0-.99D0*XI( 13-J) 

22  XP(J+12)=.01DO+.99DO*XI(J) 

C 

Q  ******************************************************************** 

C  CALCULATION  OF  THE  UNSTEADY  PRESSURE  DISTRIBUTION  AT  THE 
C  LOCATIONS  SPECIFIED  BY  THE  G-L  METHOD. 

Q  ******** ************************************************************ 

c 

222  DO  444  J=»1.24 

IF(SOC-100.)  202.201.201 
201  AVP=DCMPLX(O.ODO,O.ODO) 

AVJ=DCMPLX(O.ODO,O.ODO) 

AVK=DCMPLX(0.0D0,0.0D0) 

GO  TO  205 


-  EVALUATE  THE  CASCADE  SUMMATIONS  THAT  ARE  FUNCTIONS 

- OF  XP 

202  CALL  CCSUM  (5 , OMEGA, XP(J) ,XC ,SOC,EPS,ERR,AVP,IAVP) 

CALL  CCSUM  (6.0MEGA,XP(J),XC,S0C,EPS,ERR,AVJ,IAVJ) 
AVK=DCMPLX(0.D0,0.D0) 

IF(C.EQ.O.DO)  GO  TO  205 

CALL  CCSUM  (7 ,OMEGA,XP(J) ,XC ,SOC,EPS,ERR,AVK, IAVK) 

205  Z=ZED(OMEGB ,XP(J) ,ERR) 

CALL  U PDBC(TPD(J?,PDU(J),PDV(J), OMEGA, OMEGB,XP(J), EPS, SOC.YM.AM, BE 
1TA) 

444  CONTINUE 

********************************************************************* 
INTEGRATION  OF  UNSTEADY  PRESSURE  ALONG  THE  CHORD  TO  FIND  UNSTEADY 
LIFT  AND  PITCHING  MOMENT. 

-  BY  THE  GAUSS-LENGENDRE  QUADRATURE  METHOD 

********************************************************************* 

IF(M.EQ.l)  GO  TO  602 
CLU“DCMPLX(0 .DO ,0 .DO ) 

CLV-DCMPLX(0 . DO , 0 . DO ) 

CMU=DCMPLX(0.D0,0.D0) 

CMV*DCMPLX(0 .D0,0.D0) 

DO  26  J-1.12 
A* .01D0*W(13-J) 

CLU*A*PDU(J)+CLU 

clv-a*pdv(j)+clv 

CMU-A*PDU(J)*XP(J)+CMU 
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26  CMV«A*PDV(J)*XP(J)+CMV 
DO  28  J-13.24 
B=.01D0*W(J-12) 

CLU«B*PDU(J)+CLU 
CLV-B*PDV(J)+CLV 
CMU-B*PDU ( J ) *XP( J )+CMU 
28  CMV-B*PDV(J)*XP(J)+CHV 
M=M+1 
GO  TO  12 

602  DO  37  J-1,12 
A».99D0*Wil3-J) 

CLU-A*PDU(J)+CLU 
CLV=A*PDV(J)+CLV 
CMU“A*PDU  (  J  )  *XP(  J  )+CMU 
37  CMV-A*PDVm*XP(J)+CMV 
DO  39  J-13,24 
B= . 99D0*W( J-l 2) 

CLU=B*PDU(J)+CLU 
CLV«B*PDV(J)+CLV 
CMU=*B*PDU(J)*XP(J)+CMU 
39  CMV=B*PDV(J)*XP(J)+CMV 

******************************************************************** 
FIND  THE  COEFFICIENTS  OF  UNSTEADY  LIFT  AND  PITCHING  MOMENT 
AND  LOCATION  OF  UNSTEADY  CENTER-OF-PRESSURE 
******************************************************************** 


CLU=CLU/(2.D0*PI) 

CLV-CLV/(2.D0*PI) 

CMU— CMU/(4.D0*PI) 

CMV—  CMV/(4.D0*PI) 

CL=CLU+CLV 

CLWH=CL*EXL 

CLX*DREAL(CL) 

CLY-DIMAG(CL) 

clmag»dsqrt( clx**2-k:ly**2  ) 

PHASE=57.29578D0*ATANF(CLX,CLY) 

CLWX=DREAL(CLWH) 

CLWY-DIMAG(CLWH) 

PHAWH-57.29578D0*ATANF(CLWX,CLWY) 

CM-CMU+CMV 

CMX-DREAL(CM) 

CMY=DIMAG(CM) 

CMMAG=DSQRT(CMX**2+CMY**2) 

PHASM- 5 7 . 2 957 8D0 *ATANF ( CMX , CMY ) 

XCP- . 5D0-CMMAG/CLMAG 

WRITE(6 ,800)  CL, CLMAG , PHASE , CLU , CLV , CLWH , PHAWH , CM, CMMAG , PHASM, CMU , 
+CMV.XCP 

800  FORMAT(/5X. 'UNSTEADY  LIFT  W.R.T.  MIDCHORD' /5X, 'L"  /(2*PI*RHO*W 
1*WD*EXT)7/ 


1 12X, 'CL"' .12X. 'CLMAG' ,4X. 'PHASE (DEG) ' .1 IX, 'CLU 
l'CL*  W.R.T.  LEADING  EDGE  ,4X.' PHASE  (DEG) '/ 


,14X, 'CLV"',  5X, 


a  uu  n«nt  a  •  uuiwihu  uvvu  • -rav  (  *  / 

15X.F8.4, F8.4,2X,F8.4,4X,F8.4,5X,3(F8.4,'  ',F8.4,2X),5X,F8.4// 
15X, 'UNSTEADY  PITCHING  MOMENT  COEFF.  W.R.T.  MIDCHORD  / 

15X  '  M" /(4  •  *DT*BUn*U*Lm*ITYT  \  '  /  /  1  5Y  'r M"'  1  lY  'r MMifi'  LY  #PUACP  ^nCi 

1)'’,9X  ,'CMU" 

1,''  £8.4,2 
1/C  =’,F8.4) 


PHASE(DEG 

_.4 

XCP 


******************************************************************** 


PRINT  OUT  UNSTEADY  PRESSURE  DISTRIBUTION,  IF  JRUN  -  1 

******************************************************************** 
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IF(JRUN.EO.O)  GO  TO  25 
ANG1-57 .2957§DO*DATAN(-OMEGB) 

WRITE(6,250) 

250  F0RMAT(///4X, 'UNSTEADY  PRESSURE  DISTRIBUTION' /17X, 'TPD/(RHO*W*WD)' 
1/5X. '%  CHORD’ ,6X, 'REAL' ,6X, 'IMAG' ,6X, 'MAG' ,5X, 'PHASE(DEG)' , 
12X,TPU/(RH0*W*WD)' ,3X,  PUMAG',4X,'TPV/(RH0*W*WD)',3X,'PVMAG', 

UX, 'PHASE  WRT  L.E.'/) 

DO  75  J-1,24 
XOC=50.*(1.+XP(J)) 

CPD-TPD(J) 

TPDR-DREAL(CPD) 

TPDI-DIMAG(CPD) 

PMAG=DSQRT(TPDR**2+TPDI**2) 

ANG=*57 .29578DO*ATANF  (TPDR,  TPDI ) 

ANGL-ANG+ANG1 
PUMAG»CDABS(PDU£J j  ^ 


PVMAOCDABS  (PDV  ( 


WRITE(6 ,350)  XOC.CPD.PMAG.ANG.PDUU)  ,PUMAG,PDV( J)  .PVMAG.ANGL 
350  FORMAT(4X,F8.4.2X.2(F8.4,2X),F8.4,4X,F8.4,2X,2F8.4,lX,F8.4,lX,2F8. 


14,1X,F8.4,1X,F8.4) 
75  CONTINUE 
25  CONTINUE 
20  CONTINUE 
STOP 
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SUBROUTINE  UPDBC  (TPD.PDU, PDV, OMEGA, 0MEGB,XP, EPS, SOC.YM.AM, BETA) 
******************************************************************** 

ACCORDING  TO  EQUATION  (39)  IN  SHEN'S  THESIS 

-  CALCULATION  OF  UNSTEADY  PRESSURE  FUNCTION 

******************************************************************** 

IMPLICIT  REAL*8(A-H,0-Z,$) 

COMPLEX*!. 6  Z . RI.TMAU ,TMAV ,TMA , DEL ,MC , NC . EXW.DELU ,DELV , TMB ,TMBU,TMB 
1 V , TMCU , TMCV ,TMDD, TMDU , TMDv , SVOR , TMEU , TMEV ,P  DU , PDV , TPD ,HANK0 , HANK! , 
1EAVA.EAVB .AVL ,AVM,AVP, AVJ ,AVK 
COMPLEX* 16  ZED 

COMMON  /MAIN/  EAVA,EAVB.AVL,AVM.AVP,AVJ,AVK.Z,XC,ERR 

COMMON  /DELTA/  RJAO.RJaI .RJA2 ,RiO,RYl  ,RJB0,RJB1 ,  RJ  B2 ,  HANKO ,  HANK1 , 

1  MC.NC.DEL,DELU,DELV,RI,EXW,PHI,SINPH,COSPH 

DATA  PI/3.141592?DO? 

SXP-DSQRT(1 .0D0-XP**2) 

TXP=DSQRT ( ( 1 .ODO-XP ) / ( 1 .ODO+XP ) ) 

OMGAB-OMEGA/OMEGB 

TMAU-- 4 . *COSPH*YK*( (OMGAB- 1 . )*(2 . *XP*Z-SXP*RJB0)-RI*TXP*RJBl-2. *RI 
1  *OMEGA*Z/(OMEGB**2) ) 

TMAV—2 .  *SINPH*  (TXP*RJB0+2 .  *(OMGAB- 1 .  )*Z ) 

TMA-TMAU+TMAV 
TMBU— OMEGA*DELU*HANK0*TXP 
TMBV  =»-OMEGA*DELV*HANK0*TXP 
TMCU-(RI*OMEGA*DELU*AVP)/(PI*EXW) 

TMCV=*(RI*OMEGA*DELV*AVP)/(PI*EXW) 

TMDD-TXP* ( AV J-RI*OMEGA*AVM*EXW ) / (EXW*P I ) 

TMDU=DELU*TMDD 

TMDV-DELV*TMDD 

SVOR5*  2  .D0*TXP*(2  .D0*YM*(  1  .DO+XP  )+AM+(YM+AM)*AVK/  (2  .DO-EAVB)  ) 
TMEU=-COSPH*SVOR*CDEXP(-RI*OMEGB*XP) 

TMEaTMEU 

PDU=TMAU+TMBU+TMCU+TMDU+TMEU 

PDV=TMAV+TMBV+TMCV+TMDV 

TPD=PDU+PDV 

RETURN 

END 


FUNCTION  ZED  (OMEGB.XP.ERR) 

IMPLICIT  REAL*8(A-H,0-Z.$) 

COMPLEX*i 6  ZED.ZEDi.RI.RIK 
PSI=3. 14I5927DU-DARCOS (XP) 
ZED=DCMPLX(0.0D0,0.0D0) 

RI=DCMPLX(0 .0D0 , 1 . 00D0 ) 

10  CALL  BESJ  (OMEGB,K,RJK,0.00Q01D0,IEK) 
RIK=RI**K 

ZEDl «RIK*RJK*DSIN(K*PSI) 

ZED=ZED+ZED1 
ZED 1R=D REAL? ZEDl) 

ZED1I=DIMAG(ZED1) 

IF(DABSIZEDIR)-ERR)  40,40,30 
40  IF(DABS(ZED1I)-ERR)  20,20,30 
30  K*K+1 
GO  TO  10 
20  RETURN 
END 
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SUBROUTINE  CCSUM  (I,OMEGA,XP,XC,SOC,EPS,ERR,AVER,ISUM) 


******************************************************************** 
CALCULATION  OF  THE  UNSTEADY  CASCADE  SUMMATIONS  [TABLE  I  IN 
SHEN'S  THESIS] 

—BY  THE  AVERAGING  PROCEDURE. 

OUTPUT:  AVER.ISUM 

******************************************************************** 
IMPLICIT  REAL*8(A-H ,0— Z ,$ ) 

REAL*8  IC,NSC,JR1 , JR2,JK3 . JK4 .EAX .LAX.LBX 
COMPLEX*! 6  AVERR, AVERS .AVERT .AVER, 

+  TN,RT,GC,HC,TGC,THC,GGC,HHC,GHC,GL,HL,TGL,THL, 

ERE  RX 

COMMON  /AVIV/  iAU.SINXI.COSXI 
COMMON  /RSUM/  JK1 ,  JK2 , JK3 ,  JK4 


DATA  PI/3.1415927DO/ 
AVERR=DCMPLX(0.0D0,0.0D0) 
AVERS=DCMPLX(0 .0D0 ,0.0D0 ) 
AVERT=DCMPLX(0 .ODO, 0 .ODO) 

TXP=D SQRT ( ( 1 . ODO-XP ) / ( 1 . ODO+XP ) ) 

99  TN=DCMPLX(0.0D0,N*TAU) 
RT=CDEXP(TN) 

NSC=N*SOC 

IF(I .EQ.3)  GO  TO  30 

IF  (I.EQ.4)  GO  TO  40 

BIT=2.*NSC*SINXI 

RC-XC+BIT 

IC=2.*NSC*COSXI 

GC=DCMPLX(RC,-IC) 

HC=DCMPLX(RC,IC) 

TGC=CDSQRT( (GC+1 .0D0)/(GC-1 .ODO 
THC=CDSQRT ( (HC+1 .ODO ) / (HC- 1 .ODO 
IF! 1-2)  10,10.5 
5  IF (1-6) 50 ,^0,^0 
-  t=l 


i! 


C - EAVA 

10  EAX=TGC+THC-2 .ODO 
IF(I.EQ.2)  GO  TO  20 
AVER=EAX*RT 
GO  TO  100 

C - EAVB  1-2 

20  AVER-EAX 
GO  TO  100 

q _ —  AVL  1=3 

30  CALL  RSUM1 (OMEGA, NSC .EPS) 
AV  E  R=D  CMP  LX ( JK1 ,JK2)*RT 
GO  TO  100 

C - AVM  1=4 

40  CALL  RSUM2( OMEGA, NSC .EPS) 
AVER=DCMPLX(JK3,JK4)*RT 
GO  TO  100 

C - AVP  1=5 

50  EAX=DREAL(TGC+THC) 


LAX=DREAL(TGC*THC) 

ARG=EAX*DSQRT( 1 .0D0-XP**2 ) /( ( 1 .ODO+XP )-LAX*( 1 .ODO-XP ) ) 

RL= 1 .ODO+BlT 

GL=DCMPLX(RL,-IC) 

HL=DCMPLX( RL, IC ) 

TGL-CDSQRT ( (GL+1 .0D0)/(GL-1 .ODO) 

THL-CDSQRT ( ( HL+1 . ODO ) / ( HL- 1 . ODO ) 

E  BX  =D  REAL ( TGL+THL ) 

LBX=DREAL(TGL*THL) 

ARL=EBX*DSQRT(1 .0D0-XP**2 ) /( ( 1 . ODO+XP )-LBX*( 1 .ODO-XP ) ) 
AV  E  R=  2 , 0  DO  *RT*  (  D  AT  AN  (  (  ARG - ARL ) / ( 1 . 0  DO+ARG  *ARL ) ) ) 

GO  TO  100 


!» 


1980 


15300 


-34- 


January  5,  1980 
ICS  :mmj 


C - AVJ  1-6 

60  GGC-TGC/ (GC-XP) 
HHC-THC/ (HC-XP) 
GHC-GGC+HHC 
IF(I.EQ.7)  GO  TO  70 
AVER-GHC*RT 
GO  TO  100 

C - AVK  1-7 

70  AVER-GHC 


AVERGING  PROCEDURES  TO  FIND  THE  CONVERGED  VALUES  OF  SUMMATIONS 

100  AVER-AVER+AVERS 
AVERS-AVER 
IF(N)  120,110,110 
110  N— N 

GO  TO  99 
120  N— N 

IF(N-3)  170,140,130 
130  AVERT=AVERR/FL0AT(N-3) 

140  AVERR-AVER+AVERR 
IF(N-3)  170.170.150 
150  ERE-AVERR/FLOAT(N-2)-AVERT 
EREl-DREAL(ERE) 

ERE2=DIMAG(ERE) 

IF(DABS(ERE1)-ERR)  160,160,170 
160  IF(DaBS(ERE2)-ERR)  222,222,170 
170  IF(N.GE .50)  GO  TO  444 
N=N+1 
GO  TO  99 

222  IF(N- 10 >99,99,444 
444  ISUM=N 
RETURN 
END 
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SUBROUTINE  RSlMi  (OMEGA, NSC, GS I) 
IMPLICIT  REAL*8  (A-H,  0-Z.$) 

REAL *8  NSC, LA, JK1. JK2,JK3, JK4 
DIMENSION  X(20) ,W(20) 

COMMON  /AVIV/  TAU, SIN  XI,  COS  XI 
COMMON  /RSIM/  JK1 , JK2 , JK3 , JK4 
10  NNN-10 

X(1  >~0.  973 90652 85D0 
X(2)— 0.8650633666D0 
X(3)«- 0.6794095682D0 
X(4>— 0.4333953941D0 
X(5)— 0.1488743389D0 
X(6)»- X(5) 

X(7)— X(4) 

mpm) 

xao)— xu> 

X(ll)— l.ODO 
W(1 )»0. 06671 3443D0 
W(2)=0. 1494513491D0 
W(3)“0. 2190863625D0 
W(4)“0. 2692667 193D0 
W(5)-0. 2955242247D0 
W(6)«W(5) 

W  ( 7 ) =W ( 4 ) 

W(I0)=W(1) 

W(11)»0.0D0 
15  PI-3.1415926D0 
GS  -GS I*PI / 1 80 .  OD  0 
VM -OMEGA 

F=2 . 0D0*NSC*DSIN (GS ) 

E=2 . ODO*NSC*DCOS (GS ) 

LAMDA—1 

TA=0. ODO 

TB=0.  ODO 

DO  400  N=l, 200 

LAM  DA-LAM  DA+2 

IF  (LAMDA-1 )  51,51,52 

51  N4-NNN+1 
GO  TO  53 

52  N4-NNN 

53  SA-O.ODO 
SB -0. ODO 

DO  500  M-1.N4 
LA-X  W  HLAMDA+l  .ODO 
AL-(LA+F-1 . 0D0)**2+E**2 
AL 1 -2 . OD 0* (LA+F-1 . OD 0 ) 

AL2-AL**2 

AL3=AL**3 

AL4-AL**4 

AL5=AL**5 

BE«(LA+F)**2-1 . 0D0+E**2 
BE1-2. ODO*(LA+F) 

A-BE/AL 
B«2.0D0*E/AL 
H-A**2+B**2 
SQRTH-DSORT (DABS (H ) ) 

SH1«SQRTH*H 

SH2-SH1*H 

SH3-SH2*H 

C-2.0D0*(A+SQRTH) 

SORTC-DSQRT (DABS (C) ) 

SQ1C*SQRTC 

SQ2-SQ1*C 
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SQ3»S02*C 

G-SQRTC-2.0D0 

A1-(AL*BE1-BE*AL1)/AL2  ,  „ 

A2=2.0DQ*(AL*(AL-BE1*AL1-BE)+BE*AL1**2)/AL3 
AS— AL**2*(AL1+BE1) 

AT-AL*(AL1**2*BE1+2.0D0*AL1*BE) 

A3-6 . 0D0*  (AS-tAT-BE*AL  I  **3  )  /AL4 
Bl— 2. 0D0*E*AL1/AL2 
B2«4.0D0*E*(ALl**2-AL)/AL3 
B3=12.0D0*E*AL1*(2.0D0*AL-AL1**2)/AL4 
HI -2 . OD  0* (A*A1+B*B 1 ) 

H2-2. 0D0*(Al**2+A*A2+B*B2+B 1**2 ) 

H  3-2 . OD  0* ( 3 . OD  0*A 1 *A2+3 . OD  0*B 1 *B  24A*A34B*B  3 ) 

C 1  -2  .  OD  0*A  1-ffil /SQRTH 

C2-2. 0D0*A2-0. 5D0*H1**2/SH1+H2/SQRTH 

C3-2. 0D0*A3+O. 75D0*H1**3/SH2-1. 5U0*H1*H2/SH1+H3/SQRTH 

Gl-O. 5D0*C1/SQRTC 

G2— 0. 25DO*C1**2/SQ140. 5D0*C2/SQRTC 
G3-0.375D0*C1**3/SQ2-0.75D0*C1*C2/SQ1-K).5D0*C3/SQRTC 
IF  (M-NNN)  61,61,62 

61  A4A«AL2*(-AL+ALl**2+2. 0D0*AL1*BE1+BE) 

A4B«AL*AL1**2*(-BE1*AL1-3*BE) 

A4-2  4 • OD  0* ( A4A+A4B+B  E*AL 1**4 ) /AL5 

B4-4 8. OD 0*E* ( AL**2-3 . OD 0*AL*AL 1 **2+AL 1 **4 ) /AL5 

H4»2.0D0*(3.0D0*(A2**2+B2**2)+4.0D0*(Ai*A3+Bl*B3)+A*A4+B*34) 

C4A-Hl**4/SH3 

C4B-H1**2*H2/SH2 

C4C=H2**2/SH1 

C4D-H1*H3/SH1 

C4-2 . 0D0*A4-l .  875D0*C4A+4. 5D0*C4B-1 . 5D0*C4C-2. GD0*C4D-H14/SQRTH 

G4A=C 1**4/SQ3 

G4B=C1**2*C2/SQ2 

G4C=C2**2/SQ1 

G4D=C1*C3/SQ1 

G4— 15. 0D0/16. 0D0*G4A+2. 25D0*G4B-0.  75D0*G4C-G4D-K).  5D0*C4/SQRTC 
FA=G4*DC0S(Wd*LA) 

FB=G4*DSIN  (VM*LA) 

SA=SA+W(M)*FA 
SB=SB+W(M)*FB 
GO  TO  500 

62  3Z=G 
S1=G1 

52- G2 

53- G3 

500  CONTINUE 
TA-TA+SA 
TB-TB+SB 

GA— SZ*DSIN  (Wl)  /Wd-S  1*DC0S  (Ud) /tti**2+S2*DSIN  (ttd)  /VM**3+(S3*DCOS  (W1 
1)+TA)/WM**4 

GB-SZ*DCOS  (Wd)  /Ud-S  1*DSIN  (Wd) /Ud**2-S2*DCOS  (VM)  /VM**3+(S3*DSIN  (Wd) 
1+TB)/»1**4 
GB— GB 
29  CONTINUE 


IF  (LAM DA-30)  400,72,72 
72  YN-LAMDA+2.0D0 

IF  (OMEGA-O. 5D0)  511,76,76 
511  EN-12. 0D0/YN**4 

TESTA  AND  TESTB  ARE  MAXHHM  POSSIBLE  ERROR  BOUNDS  OBTAINED  FR(M 
ASYMPTOTIC  BEHAVIOR  OF  INTEGRAND 


TESTA-EN/DABSITA) 

TESTB -EN /DABS (TB ) 

LT-LAM  DA+2 

601  IF  (TESTA-0. 02D0)  75,75,400 

75  IF  (TESTB-0. 02D0)  76,76,400 
400  CONTINUE 

76  GA—  SZ  *DS IN  (Ud) /Wd-S l*DCOS  (Wd ) /Wd**2+S  2*DS IN  (Wd  )  /Wd**3+(S  3*DCOS  (Wd 
1)+TA)/Ud**4 

GB=SZ*DCOS  (Wi)  /Wl-S  1*DSIN (Wl)  /Wd**2-S2*DCOS  (Wd) /VM**3+(S3*DSIN  (Wd) 
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1+TB)/W1**4 
GB— GB 
JK2-OB 
JKl^JA 
603  RETURN 
END 


SUBROUTINE  RSUM2  (OMEGA, NSC, GS I) 

IMPLICIT  REAL *8  (A-H,  0-Z,$) 

REAL  *8  LA,MU.MU1,MU2,MU3,M2,M3,M4,M5,JK1,JK2,JK3,JK4,NSC 
DIMENSION  k(26),wt20) 

COMMON  /AVIV/  TAU , S INXI , COS XI 
CaiMON  /RSUM/  JK1, JK2.JK3, JK4 
10  NNN-10 

X(l)— 0.9739065285D0 
X(2)=-0. 8650633666D0 
X(3)=-0. 6794095682D0 
X(4)—0.4333953941D0 
1(5 )=-0. 1488743389D0 


X(6)=- X(5) 

7)— X(4) 

10)=-X(1 ) 

W(l )=0. 0666713443D0 
W(2)=0. 1494513491D0 
W(3)=0. 2190863625D0 
W(4  )=0. 2692667193D0 
W(5)=0. 2955242247D0 
W(6)=WC 


W  7 

W(85  .  . 

W(9)=W(2) 

W( 1 0 ) =W ( 1 ) 

15  PI-3. 1415926D0 

GSI  IS  CASCADE  STAGGER  ANGLE 
GS=GSI*PI/180.0D0 
CMEGA  IS  REDUCED  FREQUENCY 
VM  =ai  EGA 

NSC  IS  BLADE  NUMBER  TIMES  BLADE  SPACING  TO  CHORD  RATIO 
F-2 . 0D0*NSC*DSIN (GS ) 

NNN  IS  THE  NUMBER  OF  POINTS  USED  IN  THE  G-L  QUADRATURE 
E=2. ODO*NSC*DCOS (GS ) 

SA-0.0D0 
SB  =0 .  OD  0 
IF  (F)  26,26,27 
27  LS=1 
GO  TO  35 
26  LAM  DA— 1 
LS—F+3 

30  lamda=lamda+2 
RA-O.ODO 
RB-O.ODO 
DO  310  L-l.NNN 
LA-X  (L  )+LAMDA+l .  OD  0 
DEL«(LA+F)**2-1.0D0-E**2 
GAM  »2.0D0*E*(LA+F) 

MU»DEL**2-K5AM**2 
C-DEL/MU 
D -GAM /MU 
H-C**2+D**2 
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A=2.0D0*(C-tf>SQRT(DABS(H})) 

RA=RA+DSQRT (DABS  (A) )*DCOS (VM*LA)*W(L) 

RB  =RB-H>S  QRT  (DABS  (A)  )  *DS  IN  (W1*LA )  *W(L  ) 

310  CONTINUE 
SA=SA+RA 
SB-SB+RB 

IF  (LAMDA+2-LS)  30,32,32 
32  LS=LAMDA+2 
35  LA=LS 

DEL=(LA+F) **2-1 . 0D0-E**2 
DELl=2. ODO*(LA+F) 

GAM=E*DEL1 

GAM  1=2. 0D0*E 

MU=DEL**24GAM**2 

MU 1 =2 . OD  0  * (DEL*DEL 1+GAM  *GAM  1 ) 

MU2=2 .  ODO*  (DEL  1 **2+2 . OD  0*DEL)+8 . OD  0*E  **2 

MU3=12.0D0*DEL1 

C=DEL/MU 

C1=DEL1/MU-DEL*MU1A*1U**2 

C2=2.0D0/MU-2.0D0*DELl*MUlA'tU**2-DEL*MU2/MU**2+2.0D0*DEL*MUl**2/MU 

1**3 

C3=3.0D0*(MU**2*(-2.0D0*MU1-DEL1*MU2-4.0D0*DEL*DEL1)+2.0D0*MU*MU1* 
1  (DEL1*MU1+DEL*MU2  )-2. 0D0*DEL*MU1**3)  /MU**4 
D=GAM  MU 

D  l=GAli /MU-GAM  *MU  1  AlU**2 

D2=-2.0D0*GAM1*MU1/MU**2-GAM*MU2MU**2+2.0D0*GAM*MU1**2/MU**3 
D 3=  (-6 .  OD 0*E*MU2-GAM  3 )  /MU**2-»MU  1  * ( 1 2 .  OD 0*E*MU  1+6 .  QDO*GAM *MU2 )  M 

1U **3_6 . OD 0*GAM  *MU 1 **3 M U**4 
H=C**2+D**2 
H1=2.0D0*(C*C1+D*D1) 

H2=2. 0D0*(C1**2+D 1**2+C*C2+D*D2) 

H3=6. 0D0*(C1*C2+D l*D2)+2. ODO*(C*C3+D*D3) 

SQRTH=DS QRT (DABS (H ) ) 

A=2.0D0*(C+SQRTH) 

Al=2. ODO *C1+H1/SQRTH 

A2=2 . 0D0*C2-0 . 5D0*H 1**2 / (H*SQRTH )+H2 /SQRTH 

A3=2.0D0*C3-H).75D0*Hl**3/(H**2*SQRTH)-1.5DO*Hl*H2/(H*SQRTH)+H3/SQR 

1TH 

SQRTA=DSQRT (DABS (A) ) 

G=SQRTA 

G1=0. 5D0*A1/SQRTA 

G2=-0. 25D0*Al**2/(A*SQRTA)+O. 5D0*A2/SQRTA 

G3=O.375DO*Al**3/(A**2*SQRTA)-0.75D0*Al*A2/(A*SQRTA)+O.5DO*A3/SQRT 

1A 

SN=DSIN  (VM*LA) 

CN=DCOS (Wi*LA) 

TA=-G*SN/VM-Gl*CN/«<**2-KJ2*SN/VM**3-K;3*CN/»l**4 

TB=G*CN/«>1-G1*SN/W1**2-G2*CN/»1**3-K;3*SN/»1**4 

LAMDA=LS-2 

LF=LS+20 

UA=O.ODO 

UB=O.ODO 

DO  400  N =1,200 

LAi  1  DA  =L  AM  DA+2 

VA=O.ODO 

VB-O.ODO 

DO  500  M-l.NNN 

LA  =X  CM  )  +LAM  DA+1 .  CD  0 

DEL=(LA+F)**2-1. 0D0-E**2 

DEL  1=2. 0D0*(LA+F) 

GAM=E*DEL1 

GAM 1=2.  ODO*E 

MU=DEL**2+GAM**2 

M2«MU**2 

M  3*tl  2*MU 

M4»tl  3*MU 

M5*I4*MU 

MU  1  =2 .  OD 0*  (DEL*DEL  1+GAll  *GAM  l ) 
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HU  2-2 .  OD  0*  (DEL  1  **2+2 .  OD  0*DEL  )+8 .  OD  0*E**2 

MU3«12.0D0*DEL1 

C-DEL/MU 

C1-DEL1MU-DEL<MU1M2 

C2-2 .  ODO/MU+  (-2 .  OD  0*DEL  1  *MU  1-DEL^4U2 )  M  2+2 .  OD  0*DEL*1U  1**2  M  3 
C3A*1 2*  (-2 .  0D0*MU  1  -DEL  HMU2-4 .  OD  0*DEL*DEL  1 ) 

C3B  -2 .  OD  0*MU*MU  1  *  (DEL  1  *MU  140EL'IMU2  ) 

C3-3.  ODO*  (C3A4C  3B-2  •  ODO*DEL  *MU1  **3 )  Mk 
C4A— 2 . 0D0*M  3*  (MU2+4 .  OD  0*DEL  1  **2+2 .  OD 0*DEL ) 

C4B=*1 2*  (4 .  OD  0*1U1  **2+4 .  OD  0*DEL  141U 1  *MU  2+1 6 .  OD  0*DEL*DEL  1*MU  1+0EL*MU 
12**2) 

C4C  =MU*MU  1**2 *  (-4 .  OD 0*DEL  1  *MU  1-6 .  OD 0*DEL»MU2  ) 

C4=6 . OD 0* ( C4A+C 4B+C 4C+4 . OD 0*DEL*M  U 1 **4 ) M 5 
D-GAM/MU 

Dl-GAM 1  MU-GAM  *MU1  M2 

D2=  (-2 .  OD  0*GAM  1*MU  1-GAM  *t4U2 )  Ai  2+2 .  OD  0*GAM  *MU  1**2  M  3 

D3A=(-6.0D0*E1'MU2-GAM*MU3)M2 

D3B=6. 0D0*MU1*(2. 0D0*E^1U1+GAM*MU2)  Al  3 

D3=D 3A+D 3B -6 . OD  0*GAM*MU 1 **3  M  4 

D4A— 8 .  OD  0*M  3*  (E*MU3+3 .  OD  0*GAM  ) 

D4B=M2*(48. 0D0*E*MUl*MU2+8. 0D0*GAM*MUl*MU3+6.  ODO*GAM*MU2**2 ) 

D4C— 1 2.  0D0*MU*MU  1**2 * (4 .  OD 0*E*MU  1+3 .  OD 0*GAM *MU2 ) 

D4-  (D  4A-H)  4B-H)  4C+24 .  OD  0*GAM  *MU  1  **4 )  M  5 

H*C**2-H)**2 

Hl=2. 0D0*(C*C1+D*D 1 ) 

H2=2 . 0D0*(C1**2+D1**2+C*C2+D*D2) 
H3=6.0D0*(Cl*C24Ol*D2)+2.0D0*(C*C3+D*D3) 

H4=6. 0D0*(C2**2+D2**2  )+8. 0DQ*(Cl*C3+Dl*D3)+2.  GD0*(C*C4-H)*D4) 

SQRTH  =DS  QRT (H ) 

A«2.0D0*(C+SQRTH) 

Al-2. 0D0*C1+H1/SQRTH 

A2=2.0D0*C2-0.5D0*H1**2/(H*SQRTH )+H2/SQRTH 

A3=2? crf)&*C  3$A3A^  (0 .  75D0*H1**2/H-1. 5D0*H2)+H3/SQRTH 

A4=2^0D0*ci+A4A*(-l™)5D0*Hl**2/H+4. 5DO*H2)-(l.  5D0*H2**2+2.  0D0*H1* 
1H  3 ) / ( H  *S  ORTH )+H  4 / SQRTH 
SQRTA=DSQRT (DABS (A) ) 

G4A=A1**2/(A**2*SQRTA) 

G4=G4A* (-15. 0D0/16. 0D0*Al**2/A+2. 25D0*A2)- (0. 75D0*A2**2+A1*A3)/ (A* 
lSQRTA)+0. 5D0*A4/SQRTA 
VA=VA+G4*DC0S  (VM*LA)*U(I>1) 

VB  =VB4G4*DS  IN  (VM*LA  )  *W(ll ) 

500  CONTINUE 
UA=UA+VA 
UB=UB+VB 

GA=SA+TA+UA/VM**4 
GB =SB+TB+UB/WM**4 
GB=-GB 
29  CONTINUE 

IF  (LAMDA-30)  400,510,510 

510  YN  =L  All  DA+2  •  OD  0 

IF  (CAT EGA-0.  5D0)  511,76,76 

511  EN=12. 0D0/YN**4 
TESTA=DABS (EN/UA) 

TESTB=DABS(EN/UB) 

TESTA  AND  TESTB  ARE  THE  UPPER  BOUNDS  FOR  ERROR  OBTAINED  FR(M 
ASYMPTOTIC  BEHAVIOR  OF  INTEGRAND 
LT=LAMDA+2 

601  IF  (TESTA-0. 01D0)  75,75,400 

75  IF  (TESTB-O. 01D0)  76,76,400 
400  CONTINUE 

76  GA-SA+TA+UA/WM**4 
GB“SB+TB+UB/WM**4 
GB— GB 

JK3=<jA 
JK4CB 
603  RETURN 
END 
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SUBROUTINE  BESJ  (X,N,BJ,D, IER) 
IMPLICIT  REAL*8  (A-H,  0-Z,$) 
BJ-0.0D0 
IF  (N )  10,20,20 
10  IER-1 
RETURN 

20  IF  (X)  30,30,31 
30  IER-2 
RETURN 


31  IF  (X-15.0DO)  32,32.34 

32  NTEST-20. 0D0+10. 0D0*X-X**2/3 


GO  TO  36 


34  NTEST-9O.ODO+X/2.0DO 
36  IF  (N-NTEST)  40,38,38 
38  IER “4 


RETURN 


40  IER-0 
N1=N+1 
BPREV-0.0D0 

COMPUTE  STARTING  VALUE  OF  M 
IF  (X-5.0D0)  50,60,60 
50  MA-X+6.0D0 
GO  TO  70 

60  MA=1 . 4D0*X+60. ODO/X 
70  IX-X 

IlB=N+IX/4+2 

MZERO*MA 

IF  (tlA-MB)  80,90,90 
80  MZERO-MB 

SET  UPPER  LIMIT  OF  M 
90  MMAX-NTEST 

100  DO  190  M  *MZ  ERO  ,MM  AX,  3 
SET  F(M),F(M-1) 

R1 1  =1 .  OD  -2  8 
m=0.0D0 
ALPHA “0.0D0 

IF  (M-Cl/2)*2)  120,110,120 
110  JT—1 

GO  TO  130 
120  JT=1 
130  M2=M-2 

DO  160  K-1.M2 

MK=M-K 

HMK=2 . 0D0*D FLOAT  OIK)  *R1  1  /X-EM 
Rl-RU 


EM1-BMK 

IF  O-IK-N-1)  150,140,150 
140  BJ=BMK 
150  JT  =- JT 
S-l+JT 

160  ALPHA  =ALPHA+BtlK*S 
K-IK=2.0D0*RU/X-R1 
IF  (N)  180,170,180 
170  BJ-BMK 
180  ALPHA -ALPHA4BMK 
BJ-BJ /ALPHA 

IF  (DABS (BJ-BPREV)-DABS (D*BJ) )  200,200,190 
190  BPREV-BJ 
IER-3 

200  RETURN 


END 
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SUBROUTINE  BESY  (X.N.BY.IERJ 
C  CHECK  FOR  ERRORS  IN  N  AND  X 

IMPLICIT  REAL *8  (A-H,  0-Z,$> 

IF  (N)  180,10,10 
10  IER-0 

IF  (X)  190,190,20 

C  BRANCH  IF  X  LESS  THAN  OR  EQUAL  4 

20  IF  (X-4.0D0)  40,40,30 
C  CCMPUTE  YO  AND  Yl  FOR  X  GREATER  THAN  4 

30  T1-4.0D0/X 
T2*T 1*T 1 

P0-( ( ( (-0. 0000037043D0*T2+0. 00001 7356 5D0)*T2-0. 000048761 3D0)*T2+0. 
10001 7343D0)*T2-0. 001 753062D0)*T 2+0. 3989423D0 

Q0-( ( ( (0. 0000032312D0*T2-0. 0000142Q78D0)*T2+O. 0000342468D0)*T2-0.0 
1 0008697 9 IDO )*T 2+0. 0004564324DO)*T2-0.01246694DO 
Pl-((((0. 0000042414D0*T2-0. 0000200920DO)*T 2+0. 0000580 759D0)*T2-0.0 
100223203D0)*T2+O.002921826D0)*T2+O.3989423D0 

Ql»( ( ((-0.0000036594D 0*T 2+0. 0000 162 2DO)*T 2-0. 0000398708D0)*T2+0.(X) 
101064 74 lD0)#T2-0.0006390400D0)*T2+O.03740084D0 
A»2. ODO/DSQRT (X) 

B-A*T 1 

C-X-0.7853982D0 
Y0-A*P0*DSIN (C)+B*QO*DCOS (C) 

Yl— A*P1*DC0S(C)+B*Q1*DSIN  (C) 

GO  TO  90 

C  COMPUTE  YO  AND  Yl  FOR  X  LESS  THAN  OR  EQUAL  TO  4 

40  XX-X/2.0D0 
X2=XX*XX 

T”DL0G(XX)+0. 5772157D0 

SU1-0.0D0 

TEK1-T 

YO-T 

DO  70  L-1,15 

IF  (L-l )  50,60,50 
50  SUl=SUM+l . ODO/DFLOAT (L-l ) 

60  FL-L 
TS-T-SIN 

TERM - (TERM * (-X2 ) /FL**2 )*( 1 . OD 0- 1 . OD 0/ ( FL*T S ) ) 

70  YO-YO+TERM 

TERM  “XX* (T-0 . 5D  0 ) 

SIM-O.ODO 

YlwrERM 

DO  80  L-2,16 

SlW-SlM+1 .  ODO/DFLOAT  (L-l  ) 

FL°L 

FL1-FL-1.0D0 

TS-T-SIM 

TE»1 -(TERit * (-X2  ) / (FL 1  *FL ) )  * (  (TS-O.  5D0/FL ) /(TS+0. 5D0/FL 1 )  ) 

80  Y1-Y1+TEK1 

PI2=0. 6366198D0 

Y0“PI2*Y0 

Y 1 — PI  2 /X+PI 2  *Y 1 

C  CHECK  IF  ONLY  YO  OR  Yl  IS  DESIRED 

90  IF  (N-l )  100,100,130 
C  RETURN  EITHER  YO  OR  Yl  AS  REQUIRED 

100  IF  (N)  110,120,110 
110  BY-Y1 

GO  TO  170 
120  BY-YO 

GO  TO  170 

C  PERFORM  RECURRENCE  OPERATIONS  TO  FIND  YN(X) 

130  YA-YO 
YB-Y1 
K-l 

140  T-DFL0AT(2*K)/X 
YC«T*YB-YA 

IF  (DABS(YC)-1.0D36)  145,145,141 

141  IER-3 
RETURN 

145  K-K+l 
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IF  (K-N)  150,160,150 
150  YA-YB 
YB-YC 
GO  TO  140 
160  BY-YC 
170  RETURN 
180  IER-1 
RETURN 
190  IER-2 
RETURN 
END 


FUNCTION  ATANF  (X,Y) 

IMPLICIT  REAL*8  (A-H,  0-Z,$) 

ATANF  COMPUTES  THE  ANGLE  WITH  INITIAL  SIDE, 
AND  TERMINAL  SIDE,  THE  LINE  CONTAINING  THE 
ANGLES  ARE  IN  RADIANS,  GREATER  THAN  OR  EQW 
100  FORMAT  (1H  , 'INDETERMINATE  ANGLE', 5X,'X  , 

IF  (DABS (X)-0. 0Q0Q1D0)  1,2,2 

1  IF  (DABS (Y)-0. 0000 IDO)  3,4,4 

3  WRITE  (6.100)  X,Y 
ATANF “0. 0D0 
RETURN 

4  IF  (Y)  5,6.6 

5  ATANF-4.7123889D0 
RETURN 

6  ATANF-1.5707963D0 
RETURN 

2  IF  (X)  7.8.8 

7  ATANF«3.1415926D0+DATAN(Y/X) 

RETURN 

8  IF  (Y)  9.10.10 

9  ATANF“6  •  283 1 853D0-HJATAN  (Y  /X) 

RETURN 

10  ATANF “DATAN (Y /X ) 

RETURN 

END 


THE  POSITIVE  X-AXIS, 
ORIGIN  AND  POINT  (X,Y) 
X  TO  ZERO,  LESS  THAN  2* 
F10. 5, 5X, ’Y  ,F10. 5) 
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